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Abstract 

This paper considers the robust ©-stability margin problem under polynomic structured real 
parametric uncertainty. Based on the work of De Gaston and Safonov (1988), we have developed 
techniques such as, a parallel frequency sweeping strategy, different domain splitting schemes, 
which significantly reduce the computational complexity and guarantee the convergence. 
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1 Introduction 

Robustness of control systems has been one of the central issues in the control community in the 
last two decades. Most of the research efforts have been devoted to the \i framework [H [21 [121 G2] 
and the Kharitonov framework [31 [TT1 119| . One of the well studied robustness analysis problem is 
the computation of robust stability margin under polynomic structured real parametric uncertainty. 
A number of different approaches have been proposed in the Kharitonov framework aimed at the 
nonconservative computation of the robust stability margin. Among these, we recall the geometric 
programming methods [20], the algorithm based on the Routh table |18| . and the domain splitting 
approach 0CE21Q2)] based on the Zero Exclusion Condition [Sj. In general, the algorithms in [T5K20] 
is more efficient than the algorithm in [9] . The main reason is that the algorithms in [181 ED] are 
essentially based on the Routh-Hurwitz criterion and thus only finite conditions need to be evaluated, 
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while the algorithm in [9] is based on the Zero Exclusion Condition and thus a frequency sweeping 
is essential. 

Even though a frequency sweeping is a necessity, an algorithm based on the Zero Exclusion 
Condition has its particular advantage when dealing with robust D-stability problems. For example 
[3], for high order control systems, a typical specification might be as follows: The closed- loop 
polynomials should have a pair of "dominant roots" in disks of given radius e > centered at 
z\ 2 = —u ± jv, and all remaining roots having real parts less than —a with a > (See Figure [H 
where z\ € T>i, Zi £ X>2j 2? = T> \ (J £>2 U ^s)- Then, a robust Instability margin problem can be 
defined as follows: What is the maximum perturbation of plant parameters such that the roots of 
the closed loop polynomial remain robustly in T> = {z 6 C : \z — Z\\ < e}|J{z € C : \z — zi\ < 
e} U{z € C : $l(z) < a}? Since the root region T> can be defined as a union of disjoint open subsets 
with complicated boundary in the complex plane, the robust P-stability problems in general cannot 
be solved by existing results in the /i framework or the algorithms in [181 [20] which are based on the 
Routh-Hurwitz criterion. For special cases that T> is simply connected and is defined via the Nyquist 
curve of certain rational polynomials /(s) = the robust D-stability problem of p(s) may be 
reduced to the robust stability problem of polynomial p(s) = p(f(s)).(h(s)) nh where is the degree 
of polynomial h(s) and then the algorithms in [18[ [20] may be applied. However, the complexity is 
increased substantially because the coefficients of p(s) may be very complex and the degree of p(s) 
is n g times of the degree of p(s) where n g is the degree of polynomial g(s) [T7] . 
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Figure 1: Robust ©-Stability 

The advantage of an algorithm based on the Zero Exclusion Condition is that it can be applied 
to the robust P-stability problem with arbitrary complicated root region V. What we only need 
to do is to verify whether the Zero Exclusion Condition is satisfied for all boundary point of T>. 
In situations where the root region T> is complicated, applying such algorithms becomes essential. 
However, the computational complexity can be very high. In particular, the growth of computations 
is exponential with the number of parameters. It has been shown that these type of problems are 
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in general NP-hard (see [13] and the reference therein). Moreover, since the frequency search can be 
discontinuous (as shown in [1]), we usually need to evaluate the Zero Exclusion Condition for many 
boundary points of T> to come up with a reasonably accurate solution. Therefore, there is strong 
motivation to develop efficient algorithms based on the Zero Exclusion Condition to tackle the robust 
D-stability problems. 

The algorithm proposed by De Gaston and Safonov [S] is based on the Zero Exclusion Condition 
and thus can be applied to the general robust ©-stability problem. However, there exist two problems. 

First of all, it is noted that the convergence of the algorithm in [9] was concluded upon an 
impractical assumption. That is, a domain can be divided fine enough to converge to a point (see 
[9] line 40 — 53 of page 156 in the proof of the Convergence Theorem). However, to satisfy this 
assumption, the computation complexity may be unacceptably high. In this paper, we have shown 
that it is sufficient to guarantee the convergence in computing the stability margin by guaranteeing 
that the distance between critical vertices converge to 0. Therefore, it is not necessary to divide a 
sub domain so many times to collapse it to a point. In contrast, what we need is to make the critical 
vertices crunch together. Thus, the computation can be reduced greatly. We provide two splitting 
schemes which guarantee this. 

Another problem with the algorithm in [9] is its inefficiency. One main hurdle is its tedious 
frequency sweeping. Consider a family of uncertain polynomials p(s,q), q £ Q where Q is the set 
of uncertain parameters. Let k m (uj,Q) = sup{k : ^ p(ju),kQ)} where p(jco,kQ) is the value set 
associated with frequency u> and perturbation bound k. The algorithms in [9] compute k m (uj, Q) 
exactly for each frequency u and compare to find the minimum as the stability margin. To the 
best of our knowledge, all frequency sweeping techniques in the literature follow this format. In this 
paper, we investigate a smart frequency sweeping strategy. We compute k m for n r > 1 frequencies 
in parallel. Domain splitting is also performed in parallel at each iteration level. Information is 
exchanged among all subdomains to determine which subdomain for which frequency should be 
eliminated from further consideration without obtaining the exact value of k m . The stability margin 
is achieved as the minimum record of the upper bounds of all the subdomains ever generated. The 
convergence rate is much faster than that of [9]. 

The paper is organized as follows. Section 2 introduces the robust D-stability problem and 
the work of De Gaston and Safonov [9]. Section 3 discusses the Convergence Theorem of [9] and 
different domain splitting schemes. Section 4 presents our Parallel Frequency Sweeping Algorithm. 
An illustrative example is given in Section 5 and Section 6 is the conclusion. 
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2 Preliminary 



It is well known that the stability problem of an MIMO system can be reduced to the study of the 
root location of a related polynomial [31 [B]. We consider a family of polynomials p(s, q) of degree n 
whose coefficients (ii{q) are continuous functions of ^-dimensional vector of real uncertain parameters 
q, each bounded in the interval [q~ ,qf]- More formally, we define 

p(s, q) := a (q) + ai(q)s + a 2 (q)s 2 H h a n {q)s n 

where q := (qi, • • • , qi) and the hypercube Q := {q : q~ < q% < qf ' ,i = 1, • • • , £} with the nominal 
parameter q° £ Q. 

2.1 Robust P-stability Margin 

Definition 1 Let T> be an open region in the complex plane and take p(s) to be a fixed polynomial. 
Then p(s) is said to be T> -stable if all its roots lie in the region T>. 

Definition 2 A family of polynomials V = {p(-,q) ■ q £ Q} is said to be robustly V-stable if all 
roots of p(s,q) lies in D. For special case when T> is the open left half plane, V is simply said to be 
robustly stable. 

Let Q C Q. Define value setp(z,Q) C C as follows: 

p(z,Q) := {p(z,q) : q € Q}. 

Define 

kQ := {q° + k(q-q°):q(£Q}. 
We first state the Zero Exclusion Condition for uncertain polynomials. 

Theorem 1 (|8j) The polynomial p(s,q) is robustly V-stable for all q € Q if and only ifp(s,q) is 
stable for some q € Q and ^ p(z, Q) for all z £ dT>. 

Let T>\, T>2, • • • ,T>n be disjoint open subsets of the complex plane and suppose V = {p(-,q) '■ 
q £ Q} is a family of polynomials with invariant degree. For each q £ Q and i € {1, • • • , N}, let 
ni(q) denote the number of roots of p(s,q) in T>i. Finally, assume that p(s,q°) has no roots in the 
boundary of T> = T>\ (J T>2 (J • • • T>n. Then each of the root indices Ui{q) remains invariant over Q if 
and only if the Zero Exclusion Condition ^ p(z, Q) is satisfied for all points in dT>. 

Definition 3 Suppose T> is an open subset of the complex plane with boundary dT>. Then, given an 
interval I C R, a mapping <!>£> : I — ► dT> is said to be a boundary sweeping function for T> if $>x> is 
continuous and onto; i.e., $x> is continuous and for each point z G dT>, there is some 5 € / such 
that $x>(£) = z. The scalar 5 is called a generalized frequency variable for V. 
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Let k m (8,Q) := sup{/c : ^ p($>j)(5),kQ)}. The robust ©-stability margin k max is given by 
kmax = inf^g/ k m (5, Q). In general, when V = |Jj T>i where T>i, I = 1, • • • ,N are disjointed open 
subsets in the complex plane, we can define ./V boundary sweeping functions & v : I\ — ► SI);, Z = 
1, • • • ,N respectively. Then the robust ©-stability margin is given by 

k max = , min inf k m (5, Q). 

1=1, ■■■ ,iV oZLli 

2.2 Domain Splitting Algorithms 

It is noted that the analysis of robustness under polynomic structured real parametric uncertainty 
can be converted into a simpler analysis problem dealing with multilinear structured uncertainty 

PUCES]. 

Definition 4 An uncertain polynomial p(s,q) = Y17=o a i{o) sl ^ s sa ^d t° h ave a multi-linear un- 
certainty structure if each of the coefficient functions ai{q) is multi-linear. That is, if all but one 
component of the vector is fixed, then ai(q) is affine linear in the remaining component of q. More 
generally, p(s, q) is said to have a polynomic uncertainty structure if each of the coefficient functions 
ai(q) is a multi-variable polynomial in the components of q. 

In general, there exists no analytic solution for computing exactly k max . However, the following 
Mapping Theorem can be applied to obtain a lower bound for k m (5, Q) for a family of polynomials 
of multi-linear uncertainty structure. 

Theorem 2 ([21]) Suppose an uncertain polynomial p(s, q) has a multi-linear uncertainty structure. 
Then 

conv p(z, Q) = conv {p(z, q 1 ),p(z, q 2 ), ■ ■ ■ ,p(z, q 2 ')}, Vz G dV 
where conv denotes the convex hull and q 1 , ■ ■ ■ , q 2,1 denotes the 2^ vertices of the hypercube Q. 

Let ki(5,Q) := min{fc : G conv p(<E>£>(<5), kQ)}. Then by the Mapping Theorem, ki is a lower 
bound, i.e., k\ < k m . 

Definition 5 ([9]) Critical vertices are those adjacent extreme points M a , Mp o/conv p(z, kiQ) such 
that G conv {M a ,Mp}. 

Definition 6 ([9]) m(a,/3) is the number of differing coordinates of two vertices q a ,q^ that are 
mapped by p(z, .) to M a ,Mp, respectively. 

Obviously that it follows from the Mapping Theorem that k\ = k m for m(a,/3) = 0, 1. For 
m(a, 0) > 2, a vertex path is defined as follows. 
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Definition 7 ([9]) 

A vertex path is any path between critical vertices M a ,Mp consisting of m(a, (3) straight-line 
segments defined by p(z, .) as q progresses from q a to q@ along the edges of the hypercube Q. 

Define 

k u (S, Q) := inf{£; : At least one of the vertex paths of conv p(&j)(5), kQ) intercepts the origin}. 

It is shown in [9] that k u is an upper bound, i.e., k m < k u . 

In general, it is impractical to compute k max over all frequencies. The techniques developed in 
[91 [151 [16] work essentially as follows. 

Choose a range of frequency [Si , S u ] C I and grid it as 

Oj'.= di-\ , j = 1, ■ • • ,n r n c (1) 

where n r > 2, n c > 1 are integers. Apply Algorithm 1 to compute an upper bound ku and a lower 
bound ki for k m (Sj, Q) such that " 3 1 < e, j = 1, • • • ,n r n c . Then an estimate of k max can be 
defined as 

k max := min k m (Sj, Q) 

3=1,— ,n r n c 

which satisfies 

min kj < k max < max k u 
3=1,— ,n r n c j=l,— ,n r n c 

with 

max j=l,- - ,n r n c k u — n\XD.j = i ... nrTlc kj ^ 

miiij^i,...^^ kj 
Algorithm 1 [9] — Computing k m (S, Q) 

• Step 1: Determine lower bound on k m . Designate the initial uncertain parameter domain, the 
re-dimensional hypercube Q, as Qu. 

• Step 2: Determine upper bound on k m . 

• Step 3: Iterate to converge lower and upper bounds to k m . Establish an iterative procedure 
with counter r = 1, 2, 3, • • • . For each iteration perform the following operations on subdomains 
Q rp where p represents the number of subdomains left in consideration after the rth iteration. 

• Step 3—1: Increment r, i.e., r <— r + 1. 

• Step 3-2: Make orthogonal cuts midway on the longer edges of each subdomain Q rw , w = 
1, • • • ,p in order that all edge length ratios remain within a factor of 2 of each other. Designate 
these two subdomains as Q rw and Q r r w +p\. 
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• Step 3-3: Obtain ki rw := ki(S,Q rw ) and h r(w+p) := ki(8, Q r (w+ P )) via Step 1. (Note: See [9] for 
handling exceptions). 

• Step 3-4: Obtain k Urw := k u (5,Q rw ) and K r(w+p) := k u (6, Q r ( w + P )) via Step 1. (Note: See [3] 
for handling exceptions). 

• Step 3-5: Repeat Steps 3-2 to 3-4 for each w = 1, • • • ,p. 

• Step 3-6: Define k ir := mm{ki rl , h r2 , ■ ■ • , k r(2p) } and k Ur := min{k Url , k Ur2 , ■ ■ ■ , k Ur{2p) , k Ur ^} 
and e r = kur kl k ' r ■ (Note: It is shown in [9] that fcj r _ 1 < ki T < k m < k Ur < k Ur l .) 

• Step 3-7: Eliminate from further consideration all subdomains Q rw , w = 1, • • ■ ,2p, whose 
associated ki rw > k Ur . Designate the number so eliminated as u and define a new p = 1p — u. 

• Step 3-8: Repeat Steps 3-1 to 3-8 until k\ T — > k m . The stop criteria is that e r is less than a 
chosen tolerance e > 0. 

Remark 1 In the above conventional frequency sweeping algorithm, the most important mechanism 
which impacts the efficiency is the elimination of subdomains whose lower bounds are greater than 
the minimum record of the upper bounds of all subdomains of the frequency being evaluated. This 
mechanism is implemented in Step 3-7. It has been demonstrated in \1J$ that, although the the- 
oretical increase in subdomains should be exponential, in practice the growth can be linear due to 
such mechanism. We would like to note that the elimination processes for different frequencies are 
independent and hence such independent feature leaves room for a substantial reduction of the growth 
of subdomains. 

3 Domain Splitting and Convergence 

One of the most important requirement of an algorithm is on its convergence. For example, for the 
above algorithm, it is expected that given any tolerance e > 0, the above algorithm stops at finite 
iteration, i.e., r < oo for each frequency. In this section, we investigate how a domain splitting can 
affect the convergence. 

3.1 An Impractical Assumption 

The convergence of the above algorithm was addressed in [91 [10] and a Convergence Theorem was 
proposed. However, in the proof of the Convergence Theorem^, the convergence was concluded 
upon the assumption that each subdomain converges to a single point by subdivisions (see [9], lines 
40 — 53 of page 162). In another paper [T2], the convergence was also concluded by assuming that 
a subdomain is divided fine enough (see the last paragraph in page 767 of [H]). In fact, such an 
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assumption is in general impractical to be satisfied. This is because the computation is usually very 
high to divide a domain fine enough to collapse it to a point. 

In the above algorithm, the criterion adopted in splitting a domain is that "make orthogonal cut 
midway on the larger edges" . It is also addressed in [10j [16] that a splitting of a domain should be 
made in a way guaranteeing two critical vertices remained in different subdomains. In general, there 
is more than one way to satisfy these two criteria. We would like to note that, in general, a splitting 
scheme which just consists of these two criteria is not sufficient to obtain a sequence of lower bounds 
(or upper bounds) converging to k m or a sequence of subdomains converging to a single point in R . 

For example, consider a hypercube Q = {q G R 5 : qi G [0, 1], i = 1, • • • , 5}. Let Q\ = Q. Based 
on the above two criteria, Q r can be splitted as Q r +i and Q' T+ i with k m (S, Q r +i) = k m (5, Q) at the 
r-th splitting, r = 1, 2, • • • . We cannot exclude the possibility that there exists r c < oo such that the 
following are true. 

• Critical vertices differ in coordinates q±, q2, <?3, <?4 for r > r c . 

• Coordinates q%, q2 are cut in round robin order for r > r c . 

Finally, we will end up with a degenerate hypercube, with q\, q2, q$ being constants and qs, q^ 
varying within intervals, i.e., a planar "box". Because (73, q\ can vary in intervals, it is possible that 
there is a gap between the upper bound k u and the lower bound ki, i.e., 3v > such that k u — k[ > v. 

Therefore, it is important to raise the following question: 

What kind of splitting guarantees the convergence? 

3.2 Guaranteed Critical Vertices Distance Convergence 

Consider a hypercube Q = {q € R^ : qi G [q^,qt], i = 1, • ' ' A}- Define a sequence (finite or 
infinite) of domains {Q r } iteratively as follows. 

• Step 1: Let Q\ = Q. Let r = 1. 

• Step 2: If the critical vertices of domain Q r , denoted by q ar and q@ r , differ in no more than 
one coordinate then the iteration process is terminated, otherwise choose i* G U r := {i : 
Q? r ¥= ?f r } and designate either {q G Q r : < ^ < max{^ r , q^}} or {q G Q r : 
mm{q^ r ,q^ r } < q h < q -ZipLi} as Q r+V 

• Step 3: Set r = r + 1. Go to Step 2. 

In general, there are more than one way to choose i* G lA r . Let Q r = {q G R^ : qi G [g~ r , q^ r ], i = 
1, • • • ,£}, we can define a splitting scheme as follows. 
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Definition 8 A maximal-cut is a partition of Q r as above by choosing i* £ lA r such that 

qt r ~ <i7 r = max itr - <l7v 

Another splitting scheme adopted in [10] was that the cut should be made over the coordinate that 
has been subdivided the least number of times. More formally, we define a fair-cut scheme as follows. 



Definition 9 A fair-cut is a partition of Q r as above by choosing i+ G U r such that 

_+ ..- + 



mm 



a- —a- i&A r a — a- 
Now we discuss the properties of the above two domain splitting schemes. 

Theorem 3 Let {Q r } be a sequence of domains generated as above by applying the maximal-cut 
scheme in each splitting. Then, we have that either {Q r } is a finite sequence, i.e., 3ro < co such 
that the critical vertices of Q ro differ in no more than one coordinates, or {Q r } is an infinite sequence 
such that 

lim \\ P ($ v (8),q a r)-p(<!>v(5),qPr)\\ = 

r—foo 

and 

lim ki(5,Q r ) = lim k u (5,Q r ) = lim k m (S,Q r ). 

r— *oo r— >oo r— >oo 

Moreover, the same result follows if {Q r } is a sequence of domains generated as above by applying 
the fair- cut scheme in each splitting. 

Proof. We only need to consider the case that {Q r } is an infinite sequence. Decompose the 
coordinates index set I = {1, • • • , £} as X = Zf (J2oo where 

Tf = {i € T : [q~ , qf] is divided finite many times} 

and 

^oo = 6 I : [q~ , qf] is divided infinite many times}. 

Obviously, lim^oo \ \q ar — q^ r \ \ = for the case that Tf = cf). We only need to consider the case that 
If / (f>, 2oo 7^ (j>. Note that 3ri > such that qf = qf , q~ r = q^ , Mi G Tf, Mr > r\. Define 



Then 



C = min qf - q. 



min qf - q- = ( > 0, Vr > r\. 
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Note that 3r 2 > such that 

ltr-%r<C, ViG Too, Vr > r 2 . 

We claim that U r (~)Tf = 4>, Vr > max{ri,r2}. In fact, if this is not the case, then 3i* G U r (~)Tf 
such that 

because qf r — q~ r < £, Vi G W r P|Too. It follows that Q r is split as Q r +\ and Q r+1 by dividing 
interval [q^ r , g£ r ], which contradicts to g^ = g+ , g~. = Vi G If, Vr > ri. Thus, the 

claim is true and it follows that 

Therefore, lim^oo | \q ar — q@ r \\ = 0. Since p{z,q) is a continuous function of q, it follows that 
lim^oo ||p($x»((J), g" r ) -p($x>(^),? /3r )|| = 0. By the definition of ^ and k u , we have 

lim ki(S,Q r ) = lim k u (S,Q r ) = lim k m (S,Q r ). 

r— >oo r— »oo r— >oo 

Similarly, to show that the same result follows if {Q r } is a sequence of domains generated as 
above by applying the fair-cut scheme in each splitting, we only need to consider the case that 
If / <j), Too 7^ 4>- Note that 3r 3 > such that qf r = qf rz , q~ r = q~ rz , Vi G T/, Vr > r 3 . Define 



Then 



ieXf a — a- 



1i ~ 1i w ^ 

max —r — = n s < oo, Vr > r 3 . 

ieX f q + 



it, r ^>i,r 



Note that 3r 4 > such that 

4 - il 



> n s , Mi G Too, Vr > r±. 

We claim that U r (~)If = <f), Vr > max{r3,r4}. In fact, if this is not the case, then G U r f]If 
such that 

- mm — ■ < n s 



because -| — > n s , Vi G Urf]!^. It follows that Q r is split as Q r +i and by dividing 

interval [<j£ r , <z£ r ]) which contradicts to g£ r = g^, g^. = g~, 3 , Vi £ If, Vr > r 3 . Thus, the 
claim is true and it follows that 

- g^|p = ^ (g* - gf ) 2 < £ (g+ - g" Vr > max{r 3 ,r 4 }. 

2 G -2-oo i £-^"oo 
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Therefore, by the same argument as in the maximal-cut schemes, the result follows. 



□ 



Remark 2 From the proof of the theorem we can see that both domain splitting schemes guarantee 
\\q ar — q^ r \\ — > while allow Q r — ► Qoo where Qoo is not a single point in R/. Clearly, to make a 
subdomain converge to a single point requires much more computational effort than to make \\q ar — 
q^ r \\ — ► 0. As we can see later, \\q ar — q^W — > leads to the existence of a sequence of lower bounds 
(or upper bounds) converging to k m . Therefore, an algorithm based on the maximal-cut (or fair- cut) 
splitting scheme will reduce much computational effort in computing k m than other algorithms based 
on making subdomains converge to a single point in R^. From the proof, we can also see that the 
convergence will not follow if the domain splitting is made along the larger but not the largest edges 
of each subdomain. It was remarked in 110)1 that a fair-cut avoids the problem of getting into very 
narrow and long subdomains which can decrease the convergence speed. From the proof, we can see 
that it plays a role much more than affecting the speed of convergence. It is a sufficient condition to 
the existence of a sequence of lower bounds (or upper bounds) converging to k m . We would like to 
point out that the maximal-cut scheme has better worst case convergence behavior than that of the 
fair-cut scheme. 

To see the efficiency of the maximal-cut (or the fair-cut) splitting scheme, it is helpful to compare 
the image of the last subdomain resulted from the the maximal-cut (or the fair-cut) splitting scheme 
and the image of the last subdomain obtained by the finely subdivision. The situation is shown in 
Figure [2j 

The image of the last subdomain 
(based on maximal -cut or fair-cut) 



The image 
( based on the 




the last subdomain 
impractical assumption) 



Figure 2: The Image of the Last Subdomain 



4 Parallel Frequency Sweeping Algorithm 

In this section we shall investigate a new frequency sweeping structure. 
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4.1 The Main Root of Inefficiency 

To the best of our knowledge, no effort in the existing literature has been devoted to exploit a smart 
frequency sweeping strategy. Existing techniques are basically as follows: Choose and grid a range 
of frequency. Then calculate exactly k m for each gridded frequency. Finally, compare to find the 
minimum k m and return it as an estimate of k max . 

For complicated root region T>, the number of frequencies to be evaluated for k m would be 
substantial in order to obtain a reasonably good estimate for k max . Even the computation of k m for 
each frequency is very efficient, the overall complexity is still very high, because we need to evaluate 
k m for many frequencies. 

Thus for the sake of efficiency, it is natural to conceive a smart frequency sweeping strategy. 
More specifically, we would raise the following question, 

Is it possible to obtain the stability margin k max without tightly bounding k m (5j, Q) for each 
frequency 5j ? 

The following section is devoted to answering this question. 

4.2 Parallel Frequency Sweeping Algorithm 

Consider the same set of gridded frequencies Sj, j = 1, ■ • • , n r n c defined by ([1]) and relabel them as 

{6 u -5j){i-l) (* u -fr)(j-l) 
:= dH 1 , i = l,---,n r , j = l,---,n c . 

We are now in a position to present our Parallel Frequency Sweeping Algorithm as follows. 
Algorithm 2 — Parallel Frequency Sweeping Algorithm 

• Step 1: Initialize. Set j = 1. Set k = oo. Set tolerance e > 0. Set maximal iteration number 
IT. 

• Step 2: Update k and record the number of iterations r(j) for frequency 5y by the following 
steps. 

— Step 2-1: Let Uij = {Qk} = Q, * = !,-•• , n r . Set r = 1. 

— Step 2-2: If r = IT + 1 or Uij is empty for any i € {1, • • • , n r } then record r(j) = r and 
go to Step 3, else do the following for all i such that Uij is not empty. 

* Choose Q to be any element of Uij with 

h (Sij , Q) = mm h (Sij ,Q k ). 

QkSUij 

* Partition Q into Q a and Qi, by applying a maximal- cut. 

* Remove Q from Uij. 
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* Update 

k = mm{k, k u (Sij,Q a ), k u (5ij,Q b )}. 

* Add any Q € {Q a , Qb} with two or more critical vertices to Uij. 

* Remove from Uij any Q with 

k 



(2) 



^ conv p 



1 + e 



(3) 



— Step 2-3: Set r = r + 1 and go to Step 2-2. 
• Step 3: If j = n c then STOP, else set j = j + 1 and go to Step 2. 

In Algorithm 2, n r branches of frequency sweeping are performed in parallel with starting frequen- 
cies Sa, i = 1, • • • ,n r and step size ■ Each branch of frequency sweeping is not independent, 
they exchange information. The information is applied to determine the subdomains to be eliminated 
from further consideration and to update k. Finally, k is returned as the robust stability margin. 
Algorithm 2 is visualized in the following Figure [3l 
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Figure 3: A Picture of Parallel Frequency Sweeping Algorithm. N = 3, K = 4. 



Remark 3 As we can see from Step 2-2, there are two mechanisms which contribute to the efficiency 
of the Parallel Frequency Sweeping Algorithm. First, any subdomain Q that satisfies condition (fjP, 
which is equivalent to 

k l (S ij ,Q)< T ^-, (4) 

will never be partitioned again and thus can be eliminated from further consideration. Second, any 
subdomain Q with critical vertices differing in no more than one coordinates will never be partitioned 
again and thus can be eliminated from further consideration. 

We would like to note that the proposed Parallel Frequency Sweeping provides substantial improve- 
ment on efficiency than the algorithms in J2|/. This can be explained by the significant relaxation of 
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the condition for eliminating a subdomain from consideration. By p|), |5p and we can see that 
k is the minimum record of the upper bounds among all subdomains evaluated (no matter it belongs 
to the same frequency or not) and is contracted to In contrast, in algorithms of the minimum 
record of upper bounds is obtained for the frequency being considered only. Therefore, the condition 
for eliminating a subdomain from consideration is much looser than its counterpart of algorithms 
in f^. Consequently, such a significant relaxation results in a substantial decrease of the number of 
subdomains needed to be evaluated. 

Remark 4 In Algorithm 2, at each iteration, only the domain with the smallest lower bound is 
partitioned. This mechanism differs from that of Algorithm 1 in which all domains are partitioned 
and thus the number of domains increases rapidly. We can see that Algorithm 2 effectively controls 
the growth of the number of subdomains and thus is much efficient than the conventional algorithm. 

Remark 5 It is important to note that Algorithm 2 involves only one CPU processor. It is funda- 
mentally different from the parallel algorithms which involves more than one CPU processor. 

Remark 6 The speed of computation depends on the choice of integers n r and n c . When the total 
number of frequencies is fixed (i.e., n r n c is constant), the number of branches of frequency sweeping 
n r should not be too small. Small n r will hinder the improvement of efficiency benefited from the 
parallelism. However, very large n r will not result in optimal performance either. 

In addition to the novel frequency sweeping strategy, another character of Algorithm 2 is that 
there is no tolerance criteria directly forced on the final result, however, the final result falls into 
tolerance automatically. 

Theorem 4 Suppose that the maximum iteration number IT = oo. For arbitrary tolerance e > 0, 
Algorithm 2 stops with a finite number of domain splittings for each j, i.e., r(j) < oo, Vj. Moreover, 
the final k satisfies 

r, - k k max 

< — ~ < e. 

kmax 

Proof. We first show the final k > k max . Let k u be the upper bound of domain Q C Q which 
ever appeared during the execution of Algorithm 2. Let 5 be the associated frequency of Q. Note 
that € p($>d(5), k u Q) C p(<I>£>(<5), k u Q) and thus k u > k max . Note that the final k is the minimum 
record of all such fc u 's, thus k > k max . 

We next need to show that Algorithm 2 stops with a finite iteration number r(j) for each j. 
Suppose 3j such that r goes to oo. Then 3Sij such that r goes to oo. Therefore, we can construct 
a sequence of nested domains {Q r } such that Q\ D Q2 D • • • D Q r ^ Qr+i 3 " " " • Thus by 
Theorem [3] we have that Ve > 0, 3ro < 00 such that k u (5ij,Q r ) — ki(5ij,Q r ) < j^k max , Mr > tq. 
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Thus mm{k,k u (5ij,Q ro )} - h(5ij,Q ro ) < j^k max . Note that fcj(<%, Q ro +i) > h($ij,Qr ) because 
Qr +i C Qr - Also note that k never increases, thus we have k — ki(5ij,Q ro +i) < Yr-kmax- Note 



that k > k max , we have 

e - k 
k - ki(5ij,Q rQ+ i) < Y^T e k ==> T+~e < ^»"o+i)> 

which implies that ^ conv (^p ^ j-^jQro+ij) ■ Therefore, by Algorithm 2 Q ro +i will not be 

splitted. This is a contradiction. So, we have shown that Algorithm 2 stops with a finite number of 
domain splittings for each j. 

Note that 35jj such that k m (5ij, Q) = k max . Moreover, El*?* G Q such that p($>-p((Jy ), k max q*) = 0. 
Since Algorithm 2 stops with a finite number of domain splittings for each j, we have that all 
subdomains ever generated are finally eliminated from consideration. Thus, there must exists Q* 
which contains q* be eliminated from consideration at a certain level of splitting. 

Assume that k = k when Q* is eliminated from consideration. Then either £ p($>x>(5ij), j^iQ*) 
or the critical vertices of Q* differ in no more than one coordinates. If the first case is true, then 
by p($T>(Sij),k max q*) = and q* G Q*, we have that G p($x>(<%), k max Q*). Thus by £ 
p($x>(5ij), y^Q*), we have < kmax- Obviously, the final k < k and thus 

k k — k 

i "'max 

^ "'max r ~ < 6. 



1 + e 

If the latter case is true, then we have 

ki(J>ij,Q ) = k u {5{j,Q ) = k max = k. 
The proof is thus completed. 



□ 



Remark 7 From the proof, we can see that the existence of a sequence of upper bounds converging 
to kmax is due to the convergence of the distance of critical vertices. 

Remark 8 By specifying e in Algorithm 2, we can obtain an estimate k such that < k z kmax < e. 
In this sense, we can say that Algorithm 2 provides a global solution for searching k max - However, 
it should be noted that it is not necessary a global solution for the exact robust V-stability margin 
kmax- This is because it is impossible to search the whole range of the generalized frequency 5. It is 
only feasible to perform the search over a set of discrete values of 5. 

Theorem 5 Suppose that the maximum iteration number IT < oo and that Algorithm 2 stops with 
r(j) < IT, Vj. Then the final k satisfies 

q < k k max ^ 

kmax 
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Proof. Since Algorithm 2 stops with r(j) < IT, 
generated are finally eliminated from consideration, 
as for Theorem 



Vj, we can conclude that all subdomains ever 
Thus the result follows from similar argument 

□ 



5 An Illustrative Example 

Our computational experience shows that Algorithm 2 provides a significant improvement upon 
conventional algorithms for most control problems. Moreover, the improvement depends on the 
problems and can be arbitrarily good. To illustrate, we consider an example with T> chosen to be the 
open left half plane. The applications to the problems with complicated D are in exactly the same 
spirit. 

The state space equation of the linear system is as follows: 

x = A(q)x + Bu 
y = Cx 

where 



A(q) 
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with uncertain parameter q € Q = [0, 1] x [0, 3] C R 2 . We obtained a polynomial for this system as 
p(s,q) = det{sl — A{q)) which has a multilinear structure. 

The upper bound and lower bound of k m on Q is shown in Figures 4 — 7. We can see that for 
most of the frequencies the upper bounds and lower bounds are far apart and thus the importance 
of domain splitting is obvious. 

To compute k max , we uniformly grid frequency band [0.01, 15.01] and obtain 1, 500 grid frequen- 
cies as 

wj = 0.01 j, j = !,-■■ ,1500. 

In Algorithm 2, we choose the relative error e = 0.01 and n r = 30, n c = 50. The 1, 500 frequencies 
are regrouped as 



0.01 + 0.2(i - 1) +0.01(j - 1) 



1, 



,30; j = !,••■ ,50. 



We ran the program in a Sun Spark work-station. The running time is about 80 seconds. The total 
number of domains evaluated is 1,570. We obtained k = 1.4384 which is achieved at frequency 
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Fi gure 4: k m upper bound and lower bound on Q. The upper bound is plotted in dashed line and 
the lower bound is plotted in solid line. 

^20, 18 = 9.78. By Theorem HI we can concluded that 

k — k 

< - max < e = 0.01. 

~~ k 

"'max 

To compare the performance of the conventional algorithm with that of Algorithm 2, it is fair to 
choose the tolerance e r = 0.01 in Algorithm 1. We also ran the program in the same work-station. 
The running time is about 9 hours. The total number of domains evaluated is 64, 813. We obtained 
k = 1.4380 which is achieved at frequency ^9733 = 9.783. Therefore, Algorithm 2 has a speed-up of 
400 over the conventional algorithm. Moreover, the number of domains evaluated in Algorithm 2 is 
only a small fraction (which is ~ 0.0242) of that of the conventional algorithm. The number of 
domains evaluated in Algorithm 2 and the conventional one for each frequency is shown respectively 
in Figure [8] and Figure EJ 

We can see that Algorithm 2 provides much superior performance than the conventional algo- 
rithms. The improvement comes from the characteristic eliminating mechanisms in Algorithm 2. 
More formally, we describe the eliminating process in Algorithm 2 as follows. 

Let U be a record of the global upper bound achieved by frequency oj. Let Qij C Q be a domain 
associated with frequency uij. When Qij is eliminated, i.e., ki(uJij,Qij) > is satisfied, there are 
only three cases as follows. 




17 




5 10 15 

Frequency 



Fi gure 5: k m upper bound and lower bound on Q. The upper bound is plotted in dashed line and 
the lower bound is plotted in solid line. 

• Case (i): Uij < uj. We call the elimination as Backward Pruning. 

• Case (ii): ujij > uj. We call the elimination as Forward Pruning. 

• Case (hi): uj^ = uj. We call the elimination as Present Pruning. 

All the above three types of pruning processes play important roles in Algorithm 2. However, 
there is only Present Pruning in the conventional algorithm. Therefore, Algorithm 2 has a much 
powerful pruning mechanism and is much more efficient. 

In this example, we have 24 records which are shown in Figure [10J The effectiveness of the three 
types of pruning processes are shown respectively in Figures [T2l [T3l and Q31 

6 Conclusion 

We have developed techniques such as, a parallel frequency sweeping strategy, different domain 
splitting schemes, which significantly reduce the computational complexity and guarantee the con- 
vergence. Our computational experience shows that Algorithm 2 provides a substantial improvement 
on efficiency in comparison with the conventional algorithms. 
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Fi gure 6: k m upper bound and lower bound on Q. The upper bound is plotted in dashed line and 
the lower bound is plotted in solid line. 
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Figure 10: Evolution of the global upper bound. The y-coordinate represents the record value and 
the x-coordinate represents the frequency achieving it. Two consecutive records are connected by 
dashed line. 
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Figure 11: Evolution of the global upper bound. The y-coordinate represents the record value and 
the x-coordinate represents the frequency achieving it. Two consequent records are connected by 
dashed line. 
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Figure 12: Backward Pruning. The y-coordinates represents the number of domains eliminated 
the record as Case (i). The x-coordinate represents the record index. 
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Figure 13: Forward Pruning. The y-coordinates represents the number of domains eliminated by the 
record as Case (ii). The x-coordinate represents the record index. 
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Figure 14: Present Pruning. The y-coordinates represents the number of domains eliminated by the 
record as Case (iii). The x-coordinate represents the record index. 
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